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(N ' Abstract 

CL We present some new theoretical and computational results for the stationary points of bulk 

^ | systems. First we demonstrate how the potential energy surface can be partitioned into catch- 

ment basins associated with every stationary point using a combination of Newton-Raphson and 
| eigenvector-following techniques. Numerical results are presented for a 256-atom supercell repre- 

sentation of a binary Lcnnard- Jones system. We then derive analytical formulae for the number 
of stationary points as a function of both system size and the Hessian index, using a framework 
■ based upon weakly interacting subsystems. This analysis reveals a simple relation between the 

total number of stationary points, the number of local minima, and the number of transition 
states connected on average to each minimum. Finally we calculate two measures of localisation 
for the displacements corresponding to Hessian eigenvectors in samples of stationary points ob- 
tained from the Newton-Raphson-based geometry optimisation scheme. Systematic differences 
| arc found between the properties of eigenvectors corresponding to positive and negative Hessian 

"p>! ■ eigenvalues, and localised character is most pronounced for stationary points with low values of 

' the Hessian index. 

?: 

1 Introduction 

P , 

q ! Stationary points of a potential energy surface, V, correspond to structures where the gradient of 

the potential vanishes, while the Hessian index, /, is defined as the number of negative Hessian 
v |~j ■ eigenvalues. Stationary points with I > 1 are often referred to as saddles. A number of studies have 

recently focused upon the properties of saddle points for model supercooled liquids and glasses G12I2I 
A rigorous partitioning of the potential energy surface is possible in terms of the catchment basins of 
local minima [often termed 'inherent structures' for bulk materiaP^ using steepest-descent paths. 
The superposition approach to thermodynamics, where the total partition function is written as a 
sum over local minima, follows naturally from this scheme P^ 7 ^ 9 ^ 1 ^ 11 ^ 12 * Dynamical properties can 
also be calculated within this framework if minimum-to-minimum rate constants are available, using 
master equation* 13 ! 14 ! kinetic Monte Carlof^ 16 ! 17 ! 18 ! or discrete path sampling^ techniques. These 
schemes generally consider only local minima and true transition states, i.e. stationary points with 
Hessian index one (a single negative eigenvalue) in view of the Murrell-Laidler theorem,^ which 
states that if two minima are connected by an index two saddle, then there must be a lower energy 
path between them involving only true transition states. 

The potential energy surface itself is independent of temperature, atomic masses and the coor- 
dinate system. Direct analysis of V itself has revealed connections between the global properties of 
the surface and the interparticle forcesp^ as well as a theoretical basis for the empirical resultP^^ 
that the number of different local minima increases exponentially with sizeP^^lESI The distribu- 
tion of minima as a function of the potential energy can also be found from simulation dataj^^ and 
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this approach has now been used in several studies of bulk material ! 2 '^ l 29 l : ^ 1 l*4***t*l 1 1 1*1 33 

An approximately Gaussian form is expected for sufficiently large systems BBC3EQC 3 Distributions 
of transition states and barrier heights have also been reported for bulk models P 1 ^ 42 ^ 4 '^ 4 ^ 

Proposals to extend such ideas to consider a partition of configuration space in terms of saddle 
points have been suggested in the context of supercooled liquids and glassesP This approach has 
much in common with the instantaneous normal mode theory developed by Keyes and cowork- 
er | 45 | 46 E| w j iere the focus of attention is the spectrum of Hessian eigenvalues for instantaneous 
configurations. Theoretical contributions based upon the premise of dividing configuration space 
into saddles have app eared I 48 I 49 I 5U JM1 following the observation that the Hessian index of the sta- 
tionary points sampled in simulations approaches zero around the critical temperature of mode- 
coupling theoryfJ 2 ! 52 ! 53 ! 54 ! 55 ! 56 ! 57 ! 58 ! 59 ^ ! However, in contrast with the partitioning scheme based 
upon local minima, a well-defined procedure for using other stationary points in this way has yet 
to be fully developed. In particular, the most common mapping used is the minimisation of the 
square gradient, |Vy| 2 P^but this tends not to produce true stationary points for large systems, as 
we have explained in a previous contribution.^ In ^21 we show how an alternative approach based 
upon Newton-Raphson and eigenvector-following techniques can be used to achieve the desired 
partitioning. 

A further development of the present work is the analysis of stationary points for bulk models 
in terms of independent subsystems. In Jllwe first present analytic solutions for the combinatorial 
problem that determines the number of stationary points of Hessian index / in a system containing 
N atoms. These results complement recent theoretical developments of Keyes and coworkers^EH 
based upon a random energy model approach p 4 ! 65 ! ^ !^!! an( j should be useful in future studies that 
consider the dynamics and thermodynamics of model systems.^ We then consider the validity of 
the independent subsystem approximation by examining properties of the Hessian eigenvectors for 
the samples of stationary points reported in ^21 

Since saddles beyond index one usually make no significant contribution to the dynamics of 
small molecules it is important to demonstrate that conventional dynamical approaches are not 
incompatible with sampling of higher index saddles in simulations. Just as most of the volume of 
a hypersphere is associated with the surface when the dimension becomes high enoug h, 69 most of 
the configuration space in a large system must lie near the surface of the catchment basins for the 
local minimal This observation does not invalidate the partition of configuration space into the 
basins of attraction of local minima. Nevertheless it has been suggested that as the temperature 
increases the dynamics of a supercooled liquid might be interpreted in terms of transitions between 
configuration space associated with saddles, rather than between local minima! 1 ^ 52 ^ 53 ^^ Such a 
partitioning is now technically feasible, as shown in £|21 but it is not yet clear whether it offers any 
advantages over the simpler approach based on minima and true transition states. 

2 A New Partitioning Scheme 

In this section we describe a partitioning scheme that divides the potential energy surface into 
regions associated with each stationary point. Numerical results have been obtained for the same 
binary Lennard- Jones system of 256 atoms (205 A and 51 B) as we considered in earlier workp^and 
further details of the model can be found in that publication. Qualitative changes in behaviour for 
various properties have previously been reported for this system around the mode-coupling theory 
critical temperature of T c « 0.435, as mentioned above. 

In addition to the new searches described below, all the molecular dynamics (MD) and geometry 
optimisations were repeated for the 256-atom system to correct a systematic problem that affected 
the results in our previous report P3 The present results therefore supersede those in the latter 
paper, although none of the conclusions in the previous contribution are affected. As before, 
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every 1000th configuration from the data collection phase was saved and used as a starting point 
for the following geometry optimisations: (1) minimisation using Nocedal's LBFGS algorithm,^ 
(2) a transition state search using hybrid eigenvector- followingl 7 -^ 10 ! 73 ! (3) minimisation of |VT^| 2 
using Nocedal's LBFGS algorithm.^ The implementation of these algorithms and the convergence 
criteria were the same as in previous workP^ The new geometry optimisations in the present work 
are based upon a combination of Newton-Raphson-type steps and eigenvector-following. The step 
along eigendirection a in the Newton-Raphson-type search was taken 

_ ~ 2 9a m 

x a — i i (1) 

4(i + Vi+4$s/4) 



where g a and e 2 are the component of the gradient and Hessian eigenvalue corresponding to eigen- 
vector a, respectively. The sign of e 2 determines whether the step in direction a raises or lowers 
the energy, and leads to the well-known result that Newton-Raphson-type geometry optimisations 
can converge to stationary points of any Hessian indexf^EHEH j n present case this is precisely 
the behaviour that we require, and each Newton-Raphson-type optimisation therefore began with 
a maximum of twenty such steps in combination with a dynamic trust radius scheme for the max- 
imum step 

size EEHEa Some of these initial geometry optimisations were already converged to a 
root-mean-square gradient of less than 10~ 7 caa/ °AA- However, for configurations obtained from 
the higher temperature simulations the usual behaviour was for the number of negative Hessian 
eigenvalues to fall somewhat from the initial value and then oscillate without converging. This be- 
haviour probably arises because the geometry optimisation does not proceed systematically uphill 
or downhill when the sign of an eigenvalue changes sign. We therefore employed a fixed number 
(twenty) of these Newton-Raphson-type steps before switching to an eigenvector-following search 
for a stationary point of fixed Hessian index given by the number of negative eigenvalues at the 
last Newton-Raphson iteration. The step now becomes 

±29a (2) 



1^1(1 + ^1 + 402/4) 

plus for uphill and minus for downhill, independent of the sign of e 2 , and all these searches were 
tightly converged to stationary points of the required index. A dynamic trust radius step scaling 
scheme was again employed in this phase of the optimisation. 

The results of these calculations are collected in Tables Q 121 an d El As in previous work the 
number of different minima and transition states sampled and the fraction of negative eigenvalues 
decrease markedly around T c , but the system is only trapped in a single local minimum on the 
simulation time scale well below this temperatureP^SEnEH Most minimisations of |Vy| 2 converge 
to non-stationary-points (NSP's) rather than true stationary points (SP's) of V, especially at higher 
temperature. 

Statistics for the mean potential energy difference and displacement in configuration space 
between the starting point and the final geometry after each optimisation are collected in Table 
I2J While the displacements are practically the same for each class of geometry optimisation the 
energy differences are always ordered AV m i n > AV ts > AVq2 > AVnr- By this measure the 
Newton-Raphson-type scheme has clearly succeeded in finding 'closer' stationary points to the 
original configurations, although the displacement statistics do not discriminate between them. 
This scheme also provides a true partitioning of the configuration space into regions associated 
with all the stationary points, although the precise boundaries will d epend upon details of how the 
algorithms are implemented, such as the maximum step size allowed [ 75 l 7(j l In contrast, if steepest- 
descent paths are used to divide the space into catchment basins for local minima then the mapping 
is in principle unambiguous. This result follows because steepest-descent is defined in terms of a 
first-order differential equation, for which a uniqueness theorem applies.^ 
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3 Stationary Points for Independent Subsystems 

We now extend previous result s^^^Sl f or th e number of local minima and transition states to 
stationary points of any Hessian index. This analysis is based upon the assumption that a suffi- 
ciently large system of mN atoms can be divided into m effectively independent subsystems of A 
atoms each. Writing the number of stationary points of index I in an iV-atom system as n sp (iV, I) 
it follows that 

n sp (mN,0) = n sp (N,0) m . (3) 

The solution of this equation is 

n sp (A,0) = exp(aA), (4) 

where a is a constant. 

A similar argument can be given for the number of transition states Assuming that the 
rearrangements associated with the transition states can be localised to one subcell, the whole 
mA-atom system will be at a transition state when one of the subsystems is at a transition state 
and the rest are at a minimum, so 

n sp (mN, 1) =mn sp (N, l)n min (A) m -\ (5) 

with solution 

n sp (A, 1) = aN exp(aA), (6) 

where a is a another constant. 
For index two saddles we have 

n sp (mN, 2) = m n sp (N, 2) n sp (A, O)™" 1 + m{m ~ l) n sp (N, l) 2 n sp (A, 0) m ~ 2 , (7) 

where the two terms correspond to the two ways that an index two saddle for the mA-atom system 
can be generated. More generally 

n sp (mA,J)=m! £ {[ , (8) 

{n ,n\,-',nj} 1=0 Ul ' 

where nj is the number of the m subsystems that are at a stationary point of index I, and the sum 
is over the combinations of nj values that satisfy J2i n l = m an d J2i I n l = J- The ml/ Y\j =0 nj\ 
factor in the above sum is a multinomial coefficient, and the sum itself corresponds to the terms in 
the multinomial series expansion of 

[n sp (A, 0) + n sp (A, 1) + • • • + n sp (A, I max )] m (9) 

with the appropriate value of J, where I max is the maximum value of / for the A-atom system. 
This observation enables us to derive a general expression for n sp (A, /) by writing the multinomial 

as 

am, 



3 - sp (A, 0) + x 1 n sp (A, 1) + x 2 n sp (A, 2) + • • ■ + x /max n sp (A, I m3X ) . (10) 



The terms contributing to n sp (mA, J) correspond to the coefficient of x J , so 

1 d J 



n sp (mA, J) 



J! dx" 



J 

.1=0 



(11) 

x=0 
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We can now prove that the general solution to these equations is 



n sp (N,I) 



exp(aiV), 



(12) 



where the binomial coefficient 



aN\ 



T(aN + l) 



or 



I\(aN-I) T(I + l)T(aN — 7 + 1)! 



(13) 



for integer and non-integer values of a, respectively. Substituting for n sp (N,I) from equation ([12 
gives 



n sp (mN, J) 



1 d J 

1 d J 

J! / amN 
J! 1 J 



L/=o \ 



aN' 



exp(aiV) 



a;=0 



(1 + x) 



amN 

x=0 

exp(amiV) 



exp(amN) 

I amN\ 
\ J 



ex.p(amN), 



(14) 



which is consistent with the supposition of equation (|12|) . 

From equation © it is clear that the total number of stationary points of any index must obey 



<(mA0 = <W 
Therefore, n^(N), like n sp (A r , 0), is a simple exponential: 



< t (7V)=exp( 7 iV), 



(15) 



(16) 



where 7 is a constant. 



n*° t (A r ) can also be obtained from n sp (N,I) by summing over 7: 



exp(aiV ) J2 = 2<l7V exp(aiV ) = exp [(a + a In 2)N] . 



(17) 



Our result for n sp (N, I) therefore implies a relationship between the the total number of stationary 
points, the scaling exponent for the number of minima and the number of transition states connected 
on average to each minimum, which is n sp (N, l)/n sp (N, 0) = aN. These parameters are connected 
by the equation 

7 = a + aln2, (18) 

which we anticipate may play an important role in future energy landscape analysis of bulk systems. 
For example, the ratio between the number of transition states and the number of local minima 
must scale appropriately for thermodynamic and dynamic properties to exhibit proper extensive 
or intensive behaviourP^BEHEII 

If n sp (N,I) has the form given in equation (|12j) then the probability of choosing a stationary 
point of index 7 from all the stationary points of a system containing N atoms is 



P(N,I) 



-aN 



(19) 
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This distribution can be closely approximated by the Gaussian form 



(/ _ M )2 /2(J 2 

p(NJ)a ^S^< (20) 

with fi = Na/2 and a 2 = Na/4. Although this approximation does not give the correct scaling 
behaviour for P(N,1)/P(N,0) it is quite accurate in most other respects. The difference between 
the tails of the two distributions occurs because the Gaussian tends to zero at / — > ±00, whereas 
the binomial distribution vanishes at I=— 1 and aN+1, so that I max =aiV. The latter limit provides 
a physical constraint on a, namely, a < d, where d is the dimension of the system. 

The above results complete our previous analysis, and are in good agreement with the Gaussian 
distributions found empirically for small clusters.^ For the clusters the maximum of the distribu- 
tion occurs at / ~ N — 2, and therefore a ~ 2. Whether this result holds for bulk systems with 
periodic boundary conditions is not yet clear. The value a = 2 may reflect the open boundary 
conditions for the clusters, where there must always be N — 1 stable modes or the cluster will 
dissociate.^ 

Our results also complement the empirical observation that the average potential energy of 
stationary points with index / appears to scale linearly with /. A theoretical justification for 
the latter observation has been obtained by Shell et al. within a similar framework that considers 
independent subsystems.^ These authors employ a maximum term approximation to find the most 
probable saddle index, and equation (|20j) indicates that this should be a good approximation for 
large enough systems. In particular, they predict that the density of saddle points with fractional 
index, i = I/dN, should scale exponentially with size with a coefficient that depends on i, i.e. as 
exp [N0(i)]. Using our expression for n sp (N,I) it is easy to show that 



0(i) = a + -=lnf , (21) 




which varies between a and a + a In 2 (in the limit of large N). Our derivations also complement 
results based on the Morse rules™ which can provide upper and lower bounds for the number of 
stationary points of a given indexP^^ 



4 Stationary Points and Dynamics 

The above analysis is based on the assumption of independent, or weakly interacting, subsystems. 
Previous work on finite size effects for model glass formers provides support for this picture, includ- 
ing the result that 'a system of iV = 130 particles behaves basically as two non-interacting systems 
of half the size'P^EH In this idealised limit the procedure to calculate dynamical properties, such as 
the diffusion constant, D, is straightforward. For example, suppose that each subsystem consists of 
minima that are each connected to two neighbouring minima via two transition states. If the rate 
constant corresponding to each transition state is the same, and equal to k, then the expected wait- 
ing time between transitions is l/2k. If the minima are assumed to be arranged linearly in space, 
and each connected pair is separated by the same distance d, then the mean squared displacement 
after n steps is simply nd 2 for uncorrelated transitions. The average time for n steps to occur is 
n/2k, and so the one-dimensional diffusion constant is simply D = kd 2 . For a weakly interacting 
system of m such subsystems the dimension of the configuration space is multiplied by m. The 
expected waiting time in any one of these higher-dimensional minima is now l/2mk, because there 
are now 2m transition states connected to each minimum. However, the same diffusion constant 
is obtained as for a single subsystem, because the displacement must be averaged over m times as 
many atoms. 
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Now suppose that we try to find the diffusion constant using a kinetic Monte Carlo ap- 
proach^EEEEHI anc i calculations of minima and transition states.^ In this case a statistical rate 
theory may be used to estimate the transition state k using the properties of the stationary points 
or pathways. For the single subsystem described above the average waiting time per transition 
is then l/2k, and the correct diffusion constant will be obtained so long as both transition states 
are located for each minimum visited. Now consider the system containing m weakly interacting 
copies of this subsystem. In this case it is necessary to locate 2m transition states per minimum, 
and the waiting time becomes l/2mk, so that much smaller time increments occur in the resulting 
KMC simulation. However, the correct diffusion constant should still be obtained after averaging 
displacements over all the subsystems. 

It is clear from the above discussion that KMC simulations may become much less efficient as 
the system size increases if all the atoms play an active role in diffusion. At higher temperatures 
efficiency is not really the point, since standard MD calculations provide a straightforward route 
to D. Rather, the KMC approach is of interest because it provides a coarse-grained picture of the 
relevant processes. Sampling problems do indeed appear to be an issue for KMC simulations of 
glasses, as we will report elsewhere Nevertheless, the correct diffusion constant can in principle 
be obtained from such calculations by considering only local minima and true transition states, as 
long as the transitions are Markovian. This result is not in contradiction with the observation that 
progressively higher index saddles may be sampled as the temperature increases. For example, the 
mean Hessian index of instantaneous configurations or stationary points obtained by the geometry 
optimisation procedure of £jH can easily be calculated for one-dimensional double-well (Figure ^) 
or periodic potentials! 58 ^ 84 * The average index for m copies of such non-interacting systems would 
simply be m times larger. Hence the stationary points sampled by the Newton-Raphson procedure 
of j^Jwould include progressively higher index saddles as a function of temperature, up to a limiting 
value determined by m. It may also be significant that the index falls to zero only in the zero 
temperature limit, although linear extrapolation of data from the range 0.2 < kT/e < 0.5 would 
suggest a non-zero value (Figure Interestingly, as noted previously from simulations by Doliwa 
and Heuer^l and from simple models by Berthier and Garrahanj^ if the same data is represented 
on an Arrhenius plot a straight line results, arguing against a sharp transition in the properties of 
saddles near to the mode-coupling temperature. 

It is sometimes stated that the inherent structure approach to the dynamics of supercooled 
liquids is only relevant in the temperature range where there is a separation of time scales for vi- 
brational motion and transitions between minima! ^ 87 ^ 8 * However, the above discussion on the sys- 
tem size dependence of the residence times for local minima indicates that this separation depends 
on system size as well as the temperature. Goldstein, in his landmark paper on the application 
of the potential energy landscape to the study of supercooled liquids,^ actually argued for the 
separation of two atomic time scales. The first time scale corresponds to localised vibrations of an 
atom around its mean position, and the second to the time taken for a transition that involves a 
significant displacement of this atom. Goldstein realized that, whilst a specified atom is vibrating, 
many transitions could potentially occur between minima that involve rearrangements localised 
elsewhere in the system. In practice, the inherent structure approach to dynamics will be most 
useful when the times between interminimum transitions are long compared to the time scale for 
vibrational motion. For this reason it has been argued that the system size needs to be small 
enough for this time scale separation to be present [ 9U I 91 I 88 I 9 H 

The previous discussion also shows that analysis of dynamics using a configuration space par- 
titioned into local minima and considering only true transition states is not incompatible with 
the system 'sampling' higher index saddles. However, all the arguments are based on the notion 
of independent subsystems, and we therefore examine how well this picture might apply to the 
stationary points obtained for the binary Lennard- Jones system discussed in fusing the Newton- 
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Raphson-type procedure. Two distinct indices have been calculated for each of these points, as 
described below. 

The first 'localisation' index, L, was defined as 



37V 

L a f3 = E | C 7 C 7 
7=1 



(22) 



where is component 7 of normalised Hessian eigenvector a. For every stationary point with 
Hessian index two or higher we calculated the average value of L over the /(/ — l)/2 pairs of 
eigenvectors with negative eigenvalues, L m ; nus . The statistics for stationary points of the same index 
were found to be very similar for samples obtained from the MD runs at different temperatures, 
and so averages over all the runs are presented in Figure [21 For eigenvectors corresponding to 
motion in different regions of space (or different atoms) we would expect L = 0, while for motion 
localised on the same atoms L should approach unity. The corresponding averages over all pairs 
of eigenvectors corresponding to positive eigenvalues, L p i us , were also calculated for comparison. 
Figure |21 reveals that L p i us is around 0.55 practically independent of the index. In contrast, i m inus 
is systematically smaller, particularly when the index is less than about ten. 
We have also calculated 

^=fE( c 7) 2 ) /E(S Q ) 4 > (23) 

\ 7 / 7 

for all eigenvectors with non-zero Hessian eigenvalues. N is proportional to the participation ratioj^ 
and is expected to vary between one for localised modes to about N for delocalised modes. Figure 
121 shows the results for modes with positive and negative eigenvalues, iV p i us and N m i nus , separately 
for comparison, and again we display averages over stationary points from the different MD runs. 
Here the difference between N p \ us and iV m i nus is even more marked than for -L p i us and L m inus- 

The statistics for both L and N both indicate that the characteristic displacements associated 
with Hessian eigenvectors that have negative eigenvalues are more localised and spatially inde- 
pendent than for eigenvectors associated with positive eigenvalues. For both measures it is the 
stationary points with the fewest negative eigenvalues for which this character is most pronounced. 
These results appear to agree very well with the analysis of Shell et al, who conclude that low-index 
saddles may often be described in terms of combinations of transition states of the subsystems.^ 



5 Conclusions 

We have previously pointed out the need for rate constants and related dynamical properties, 
such as the diffusion constant, to scale correctly with system sizeP^l The total energy and its 
fluctuations are extensive quantities, while the barrier heights between local minima and true 
transition states are intensive. If 'activated' processes are defined by comparing such extensive and 
intensive quantities then one would be forced to conclude that no 'activated' processes exist at any 
infinitesimal temperature in a bulk system ElSSI We have previously argued that such a definition 
is inappropriate pSl since rate constants for well-defined geometrical rearrangements are intensive 
quantities, as are the expressions used to calculate them in standard unimolecular rate theory.^ 
In this sense all transitions between the catchment basins of local minima are 'activated', since a 
potential energy barrier is involved, although a more useful definition should probably consider the 
magnitude of kT or the available energy per degree of freedom. 

Our previous results j ^ 4 ^ 25 ! 44 * now supported by independent calculations,^ indicate that the 
relevant potential energy barriers for diffusion are generally not small compared to kT in the 
supercooled region. In fact, statistical rate theories have been successfully applied to diffusion 
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in solids for nearly fifty years, and the standard approach is based on transition state theory 
for true transition states.^ At first sight, these observations might appear to be incompatible 
with the notion that the system samples mostly higher index saddles above some temperature 
threshold. However, the analysis of ^indicates that if the potential energy surface is partitioned 
into catchment basins of local minima in the usual way then the dynamics can indeed be treated by 
considering the lowest barriers between them, which are those mediated by true transition states 
Of course, the usual caveats apply to such an analysis, namely that the transitions are assumed to 
be Markovian, and a statistical theory is usually employed to calculate the required rate constants, 
often involving a harmonic approximation. 

We have shown how an alternative partitioning scheme based upon Newton-Raphson and 
eigenvector-following geometry optimisation can successfully divide the potential energy surface into 
catchment basins associated with all the stationary points. Thermodynamic^ and dynamic F 2 ^^ 51 * 
schemes might be constructed on this basis along similar lines to the methods used for the con- 
ventional partitioning in terms of local minima. However, the latter division is simpler, and the 
theoretical tools that use it are comparatively well developed P^l 

For weakly interacting subsystems we have now solved the combinatorial problem that defines 
the number of stationary points of any given index. This analysis also reveals a simple relation 
between the parameters that determine the total number of stationary points, the number of lo- 
cal minima, and the number of transition states connected on average to each minimum. We 
have further investigated the displacements corresponding to the Hessian eigenvectors of all the 
stationary points located using the Newton-Raphson-based scheme of £jHto see whether they are 
compatible with the above analysis. The results indicate that eigenvectors corresponding to neg- 
ative eigenvalues involve displacements where fewer atoms participate than for eigenvectors with 
positive eigenvalues, and that the displacements corresponding to different eigenvectors are more 
independent. Both these characteristics are most pronounced for stationary points with low val- 
ues of the Hessian index, in agreement with previous work in which such points are considered as 
combinations of true transition states for subsystems.^ 
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Table 1: Mean total energy, E, potential energy, V, kinetic energy, KE, and kinetic equipartition 
temperature, T, for seven MD simulations of a BLJ system with number density 1.2 crT? P^l A 
cutoff of 2.5 a aa was employed in these calculations together with a shifting scheme that ensured 
continuity of the energy and gradient.^ The ± values represent one standard deviation. #min, #ts, 
#G2 and #NR are the number of distinct minima, transition states, stationary points of | VU| 2 and 
stationary points of V, found for 10 3 searches, using L-BFGS minimisation EED hybrid eigenvector- 
following transition state searchingpl minimisation of |VU| 2 using the L-BFGS algorithm, and 
Newton-Raphson-type steps, respectively (excluding permutational isomers). %SP and %NSP are 
the percentage of quenches on the |VU| 2 surface that converged to stationary points and non- 
stationary points of V, respectively (out of 10 3 total) 



run 


E 


V 


KE 


T 


#min 


#ts 


#G2 


#NR 


%SP 


%NSP 


1 


-1648.904 ±0.003 


-1723 ±3 


75 ±3 


0.196 ±0.007 


1 


188 


85 


87 


17.8 


82.2 


2 


-1599.986 ±0.005 


-1700 ±4 


100 ±4 


0.262 ± 0.005 


14 


198 


419 


379 


12.6 


87.4 


3 


-1499.959 ±0.008 


-1652 ±6 


152 ±6 


0.399 ±0.015 


592 


874 


998 


995 


2.3 


97.7 


4 


-1399.959 ±0.012 


-1602 ±7 


202 ±7 


0.530 ± 0.020 


555 


892 


1000 


1000 


1.2 


98.8 


5 


-1299.946 ±0.015 


-1536 ±9 


236 ±9 


0.619 ±0.024 


998 


1000 


1000 


1000 


1.5 


98.5 


6 


-1199.951 ±0.020 


-1481 ± 11 


281 ± 11 


0.737 ± 0.028 


1000 


1000 


1000 


1000 


2.7 


97.3 


7 


-1099.952 ±0.025 


-1427 ± 12 


327 ± 12 


0.859 ±0.032 


1000 


1000 


1000 


1000 


2.3 


97.7 
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Table 2: Mean potential energy differences, Ay, and displacements, AD, between the starting point 
and the converged geometry after searching for minima (min), transition states (ts), minimising 
|Vy| 2 (G2), and performing Newton-Raphson-type (NR) optimisation!^ The ± values represent 
one standard deviation. 



run 




A-D mm 




A As 


AVg 2 


A-DG2 


AVnr 


AD N r 


1 


49 ±2 


32 ±0 


44 ±3 


33 ±1 


48 ±2 


32 ±0 


47 ±5 


32 ±1 


2 


98 ±4 


33 ±1 


96 ± 5 


34 ±1 


96 ±4 


34 ±1 


91 ± 12 


34 ±2 


3 


147 ±5 


37 ± 1 


146 ±6 


37 ±1 


142 ±6 


37 ± 1 


121 ±21 


38 ±2 


4 


197 ±8 


56 ±1 


195 ±8 


56 ±1 


187 ±9 


56 ±1 


159 ± 22 


57 ± 1 


5 


235 ±9 


133 ±9 


232 ±9 


133 ±9 


205 ± 10 


133 ±9 


170 ± 18 


133 ±8 


6 


283 ±9 


206 ± 11 


280 ± 11 


206 ± 11 


234 ± 12 


206 ± 11 


192 ± 19 


206 ± 11 


7 


333 ± 13 


293 ± 14 


329 ± 13 


293 ± 14 


263 ± 14 


293 ± 14 


216 ±21 


293 ± 14 



Table 3: isp an d ^nsp are the fractions of negative Hessian eigenvalues found after minimising | VV| , 
split into stationary points and non-stationary points of V, respectively.^ ^nr is the corresponding 
fraction for Newton-Raphson-type searches!^ The ± values represent one standard deviation. 



run 


Z N R X 10 3 


i N SP x 10 3 


isp x 10 3 


1 


0.4 ±0.7 


0.4 ±0.7 


0.0 ±0.0 


2 


1.4 ±1.1 


1.0 ±1.0 


0.6 ±0.8 


3 


4.1 ±2.1 


3.6 ± 1.8 


3.2 ±1.8 


4 


4.7 ±2.1 


3.5 ± 1.9 


3.4 ± 1.5 


5 


13.6 ±3.8 


10.9 ±3.2 


16.4 ±3.7 


6 


21.7 ±4.3 


17.7 ±3.8 


19.5 ±4.6 


7 


29.1 ±4.8 


23.9 ±4.4 


25.8 ±3.4 
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Figure Captions 



1. Mean Hessian index as a function of temperature (canonical ensemble) in reduced units of 
kT/e for the double well potential 8e(2x 4 — x 2 ), which has a barrier height of e. The Hessian 
index calculated for instantaneous configurations and after the Newton-Raphson procedure 
of <J21give indistinguishable results. The inset shows the same data presented in an Arrhenius 
plot. 

2. Variation of L and N, parameters designed to provide insight into the localisation of displace- 
ments along Hessian eigenvectors, as a function of the Hessian index. Both parameters were 
averaged over all the relevant stationary points obtained by Newton-Raphson-type optimisa- 
tion from configurations sampled in the seven MD runs described in Tabled For L separate 
averages were calculated over all pairs of eigenvectors corresponding to negative eigenvalues 
and all pairs corresponding to positive eigenvalues within each stationary point. For N sepa- 
rate averages were compiled for Hessian eigenvectors corresponding to positive and negative 
eigenvalues, as indicated. 
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kT/e 
Figure 1: 
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